Climate Water Loss Experiment - Treatment Analysis

Savannah Weaver

2022

Packages

`%nin%` = Negate(`%in%`)
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse") # workflow and plots
if (!require("lme4")) install.packages("lme4")
library("lme4") # for LMMs
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # for p-values
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("UsingR")) install.packages("UsingR")
library("UsingR") # simple.eda model assumption checker
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # lmer model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("car")) install.packages("car")
library("car") # VIFs
if (!require("AICcmodavg")) install.packages("AICcmodavg")
library("AICcmodavg") # model selection
if (!require("MuMIn")) install.packages("MuMIn")
library("MuMIn") # model selection
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # color
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown format

Background and Goals

This data was collected June - August by Master’s student Savannah Weaver, advisor Dr. Emily Taylor, and research assistants Tess McIntyre and Taylor Van Rossum. Adult male Sceloporus occidentalis were caught across the Cal Poly campus then acclimated to 4 different climate treatments. This R file analyzes the effect of experimental climate treatments on lizard body condition, osmotic balance (plasma osmolality and hematocrit), and cutaneous evaporative water loss (CEWL).

Data

Load

Read-in data that was compiled, formatted, and checked for completeness in ‘wrangling_general’. See that file for information related to the variables.

dat <- read_rds("./data/analysis_data_experiment.RDS") %>%
  # add VPD values by tmt and trial group
  left_join(read_rds("./data/HOBO_tmt_trial_diffs.RDS"), 
            by = c('tmt', 'trial_number'='trial')) %>%
  # add the per tmt group calculation of VPD from analysis_HOBO
  mutate(VPD_kPa = case_when(substr(tmt, 1, 6) == "Hot Dr" ~ 3.82,
                             substr(tmt, 1, 6) == "Hot Hu" ~ 1.07,
                             substr(tmt, 1, 6) == "Cool D" ~ 2.50,
                             substr(tmt, 1, 6) == "Cool H" ~ 0.64))
summary(dat)
##  measurement_date        type     individual_ID     mass_g     
##  Min.   :2021-06-16   exp  :804   201    :  7   Min.   : 7.00  
##  1st Qu.:2021-07-01   rehab:132   202    :  7   1st Qu.: 9.50  
##  Median :2021-07-25               203    :  7   Median :10.60  
##  Mean   :2021-07-22               204    :  7   Mean   :10.64  
##  3rd Qu.:2021-08-14               205    :  7   3rd Qu.:11.70  
##  Max.   :2021-09-01               206    :  7   Max.   :17.40  
##                                   (Other):894                  
##  hematocrit_percent trial_number temp_tmt   humidity_tmt     SVL_mm     
##  Min.   :13.00      1:175        Hot :467   Humid:468    Min.   :60.00  
##  1st Qu.:26.00      2:203        Cool:469   Dry  :468    1st Qu.:66.00  
##  Median :32.00      3:231                                Median :67.00  
##  Mean   :31.99      4:189                                Mean   :67.74  
##  3rd Qu.:37.00      5:138                                3rd Qu.:70.00  
##  Max.   :52.00                                           Max.   :77.00  
##  NA's   :408                                                            
##                    tmt          day_n        day_factor osmolality_mmol_kg_mean
##  Cool Humid (0.6 kPa):238   Min.   : 0.000   0 :134     Min.   :295.3          
##  Hot Humid (1.1 kPa) :230   1st Qu.: 4.000   4 :134     1st Qu.:336.1          
##  Cool Dry (2.5 kPa)  :231   Median : 6.000   5 :134     Median :351.3          
##  Hot Dry (3.8 kPa)   :237   Mean   : 5.705   6 :134     Mean   :354.3          
##                             3rd Qu.: 8.000   7 :134     3rd Qu.:370.0          
##                             Max.   :10.000   8 :134     Max.   :471.5          
##                                              10:132     NA's   :414            
##  CEWL_g_m2h_mean   msmt_temp_C    msmt_RH_percent cloacal_temp_C 
##  Min.   : 7.152   Min.   :24.80   Min.   :25.52   Min.   :23.00  
##  1st Qu.:19.755   1st Qu.:26.20   1st Qu.:46.11   1st Qu.:25.00  
##  Median :24.152   Median :26.74   Median :47.88   Median :26.00  
##  Mean   :24.767   Mean   :26.72   Mean   :46.74   Mean   :25.92  
##  3rd Qu.:28.505   3rd Qu.:27.11   3rd Qu.:50.50   3rd Qu.:27.00  
##  Max.   :56.066   Max.   :29.20   Max.   :56.16   Max.   :30.00  
##  NA's   :669      NA's   :668     NA's   :668     NA's   :668    
##   msmt_temp_K      e_s_kPa_m       e_a_kPa_m       msmt_VPD_kPa  
##  Min.   :297.9   Min.   :3.219   Min.   :0.9894   Min.   :1.486  
##  1st Qu.:299.4   1st Qu.:3.504   1st Qu.:1.6464   1st Qu.:1.767  
##  Median :299.9   Median :3.620   Median :1.7411   Median :1.853  
##  Mean   :299.9   Mean   :3.620   Mean   :1.6833   Mean   :1.937  
##  3rd Qu.:300.3   3rd Qu.:3.701   3rd Qu.:1.7992   3rd Qu.:2.012  
##  Max.   :302.4   Max.   :4.194   Max.   :1.9326   Max.   :3.021  
##  NA's   :668     NA's   :668     NA's   :668      NA's   :668    
##       SMI         temp_mean_tmttrial temp_SD_tmttrial humidity_mean_tmttrial
##  Min.   : 6.747   Min.   :23.30      Min.   :0.5966   Min.   :13.75         
##  1st Qu.: 9.714   1st Qu.:24.05      1st Qu.:0.7828   1st Qu.:29.21         
##  Median :10.594   Median :24.88      Median :1.0461   Median :45.24         
##  Mean   :10.599   Mean   :29.60      Mean   :1.1513   Mean   :52.94         
##  3rd Qu.:11.390   3rd Qu.:35.05      3rd Qu.:1.5191   3rd Qu.:82.84         
##  Max.   :15.063   Max.   :36.00      Max.   :1.8447   Max.   :93.15         
##                                                                             
##  humidity_SD_tmttrial    e_s_kPa      VPD_kPa_tmttrial    VPD_kPa     
##  Min.   : 4.370       Min.   :2.859   Min.   :0.1958   Min.   :0.640  
##  1st Qu.: 6.234       1st Qu.:2.992   1st Qu.:0.7925   1st Qu.:0.640  
##  Median : 7.382       Median :3.142   Median :2.0310   Median :1.785  
##  Mean   : 8.765       Mean   :4.330   Mean   :1.9985   Mean   :2.010  
##  3rd Qu.:12.297       3rd Qu.:5.639   3rd Qu.:3.1520   3rd Qu.:3.820  
##  Max.   :19.846       Max.   :5.944   Max.   :4.0640   Max.   :3.820  
## 

Split

Make sub-dataframes without recovery data / with only recovery-related data:

dat_no_rehab <- dat %>%
  dplyr::filter(day_n %in% c(seq(0,8))) 

recovery_values <- dat %>%
  dplyr::filter(day_n == 10) %>%
  dplyr::select(individual_ID,
                end_hct = hematocrit_percent,
                end_osml = osmolality_mmol_kg_mean,
                end_SMI = SMI)

recovery_v_post_exp <- dat %>%
  dplyr::filter(day_n == 8) %>%
  left_join(recovery_values, by = 'individual_ID') %>%
  mutate(delta_osml_10_8 = end_osml - osmolality_mmol_kg_mean,
         delta_hct_10_8 = end_hct - hematocrit_percent,
         delta_SMI_10_8 = end_SMI - SMI)

recovery_v_pre_exp <- dat %>%
  dplyr::filter(day_n == 0) %>%
  left_join(recovery_values, by = 'individual_ID') %>%
  mutate(delta_osml_10_0 = end_osml - osmolality_mmol_kg_mean,
         delta_hct_10_0 = end_hct - hematocrit_percent,
         delta_SMI_10_0 = end_SMI - SMI)

Check

Dates:

unique(dat$measurement_date)
##  [1] "2021-06-16" "2021-06-20" "2021-06-21" "2021-06-22" "2021-06-23"
##  [6] "2021-06-24" "2021-06-26" "2021-06-30" "2021-07-01" "2021-07-02"
## [11] "2021-07-03" "2021-07-04" "2021-07-06" "2021-07-20" "2021-07-24"
## [16] "2021-07-25" "2021-07-26" "2021-07-27" "2021-07-28" "2021-07-30"
## [21] "2021-08-08" "2021-08-12" "2021-08-13" "2021-08-14" "2021-08-15"
## [26] "2021-08-16" "2021-08-18" "2021-08-22" "2021-08-26" "2021-08-27"
## [31] "2021-08-28" "2021-08-29" "2021-08-30" "2021-09-01"

Number of measurements for each lizard:

dat_no_rehab %>%
  group_by(individual_ID) %>%
  summarise(n = n()) %>%
  arrange(n)
## # A tibble: 134 × 2
##    individual_ID     n
##    <fct>         <int>
##  1 201               6
##  2 202               6
##  3 203               6
##  4 204               6
##  5 205               6
##  6 206               6
##  7 207               6
##  8 208               6
##  9 209               6
## 10 210               6
## # … with 124 more rows

Every lizard has 6 experimental measurements: pre-tmt, mid-tmt, post-tmt, and mass checks on each of the 3 days between mid and post-tmt.

Did any of the treatment groups inherently start out with large differences in response variables?

dat %>%
  dplyr::filter(day_n == 0) %>%
  group_by(tmt) %>%
  summarise(mean(mass_g),
            sd(mass_g),
            mean(SMI),
            mean(hematocrit_percent),
            mean(osmolality_mmol_kg_mean),
            mean(CEWL_g_m2h_mean))
## # A tibble: 4 × 7
##   tmt                  `mean(mass_g)` sd(mass_…¹ mean(…² mean(…³ mean(…⁴ mean(…⁵
##   <fct>                         <dbl>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 Cool Humid (0.6 kPa)           11.6       1.35    11.7    39.6    351.    20.9
## 2 Hot Humid (1.1 kPa)            11.6       1.75    11.5    37.9    347.    21.4
## 3 Cool Dry (2.5 kPa)             11.8       1.61    11.8    39.3    346.    20.0
## 4 Hot Dry (3.8 kPa)              12.0       1.68    11.8    38.9    347.    20.9
## # … with abbreviated variable names ¹​`sd(mass_g)`, ²​`mean(SMI)`,
## #   ³​`mean(hematocrit_percent)`, ⁴​`mean(osmolality_mmol_kg_mean)`,
## #   ⁵​`mean(CEWL_g_m2h_mean)`

There are slight differences, but overall the starting values across groups are more or less the same.

Temp & RH during (all, before and after exp) CEWL measurements:

summary(dat_no_rehab)
##  measurement_date        type     individual_ID     mass_g     
##  Min.   :2021-06-16   exp  :804   201    :  6   Min.   : 7.00  
##  1st Qu.:2021-06-30   rehab:  0   202    :  6   1st Qu.: 9.50  
##  Median :2021-07-25               203    :  6   Median :10.60  
##  Mean   :2021-07-22               204    :  6   Mean   :10.65  
##  3rd Qu.:2021-08-13               205    :  6   3rd Qu.:11.60  
##  Max.   :2021-08-30               206    :  6   Max.   :17.40  
##                                   (Other):768                  
##  hematocrit_percent trial_number temp_tmt   humidity_tmt     SVL_mm     
##  Min.   :13.00      1:150        Hot :402   Humid:402    Min.   :60.00  
##  1st Qu.:28.25      2:174        Cool:402   Dry  :402    1st Qu.:66.00  
##  Median :33.00      3:198                                Median :67.00  
##  Mean   :33.75      4:162                                Mean   :67.73  
##  3rd Qu.:39.00      5:120                                3rd Qu.:70.00  
##  Max.   :52.00                                           Max.   :77.00  
##  NA's   :406                                                            
##                    tmt          day_n     day_factor osmolality_mmol_kg_mean
##  Cool Humid (0.6 kPa):204   Min.   :0.0   0 :134     Min.   :295.3          
##  Hot Humid (1.1 kPa) :198   1st Qu.:4.0   4 :134     1st Qu.:334.7          
##  Cool Dry (2.5 kPa)  :198   Median :5.5   5 :134     Median :348.3          
##  Hot Dry (3.8 kPa)   :204   Mean   :5.0   6 :134     Mean   :352.3          
##                             3rd Qu.:7.0   7 :134     3rd Qu.:367.4          
##                             Max.   :8.0   8 :134     Max.   :445.5          
##                                           10:  0     NA's   :413            
##  CEWL_g_m2h_mean   msmt_temp_C    msmt_RH_percent cloacal_temp_C 
##  Min.   : 7.152   Min.   :24.80   Min.   :25.52   Min.   :23.00  
##  1st Qu.:19.755   1st Qu.:26.20   1st Qu.:46.11   1st Qu.:25.00  
##  Median :24.152   Median :26.74   Median :47.88   Median :26.00  
##  Mean   :24.767   Mean   :26.72   Mean   :46.74   Mean   :25.92  
##  3rd Qu.:28.505   3rd Qu.:27.11   3rd Qu.:50.50   3rd Qu.:27.00  
##  Max.   :56.066   Max.   :29.20   Max.   :56.16   Max.   :30.00  
##  NA's   :537      NA's   :536     NA's   :536     NA's   :536    
##   msmt_temp_K      e_s_kPa_m       e_a_kPa_m       msmt_VPD_kPa  
##  Min.   :297.9   Min.   :3.219   Min.   :0.9894   Min.   :1.486  
##  1st Qu.:299.4   1st Qu.:3.504   1st Qu.:1.6464   1st Qu.:1.767  
##  Median :299.9   Median :3.620   Median :1.7411   Median :1.853  
##  Mean   :299.9   Mean   :3.620   Mean   :1.6833   Mean   :1.937  
##  3rd Qu.:300.3   3rd Qu.:3.701   3rd Qu.:1.7992   3rd Qu.:2.012  
##  Max.   :302.4   Max.   :4.194   Max.   :1.9326   Max.   :3.021  
##  NA's   :536     NA's   :536     NA's   :536      NA's   :536    
##       SMI         temp_mean_tmttrial temp_SD_tmttrial humidity_mean_tmttrial
##  Min.   : 7.317   Min.   :23.30      Min.   :0.5966   Min.   :13.75         
##  1st Qu.: 9.748   1st Qu.:24.05      1st Qu.:0.7828   1st Qu.:29.21         
##  Median :10.624   Median :29.74      Median :1.0461   Median :45.24         
##  Mean   :10.607   Mean   :29.61      Mean   :1.1502   Mean   :52.95         
##  3rd Qu.:11.348   3rd Qu.:35.05      3rd Qu.:1.5191   3rd Qu.:82.84         
##  Max.   :14.263   Max.   :36.00      Max.   :1.8447   Max.   :93.15         
##                                                                             
##  humidity_SD_tmttrial    e_s_kPa      VPD_kPa_tmttrial    VPD_kPa     
##  Min.   : 4.370       Min.   :2.859   Min.   :0.1958   Min.   :0.640  
##  1st Qu.: 6.234       1st Qu.:2.992   1st Qu.:0.7925   1st Qu.:0.640  
##  Median : 7.382       Median :4.323   Median :2.0310   Median :1.785  
##  Mean   : 8.758       Mean   :4.333   Mean   :1.9993   Mean   :2.011  
##  3rd Qu.:12.297       3rd Qu.:5.639   3rd Qu.:3.1520   3rd Qu.:3.820  
##  Max.   :19.846       Max.   :5.944   Max.   :4.0640   Max.   :3.820  
## 
dat_no_rehab %>%
  group_by(type) %>%
  summarise(mean(msmt_temp_C, na.rm = T),
            sd(msmt_temp_C, na.rm = T),
            mean(msmt_RH_percent, na.rm = T),
            sd(msmt_RH_percent, na.rm = T),
            mean(msmt_VPD_kPa, na.rm = T),
            mean(msmt_VPD_kPa, na.rm = T))
## # A tibble: 1 × 6
##   type  `mean(msmt_temp_C, na.rm = T)` sd(msmt_temp_C,…¹ mean(…² sd(ms…³ mean(…⁴
##   <fct>                          <dbl>             <dbl>   <dbl>   <dbl>   <dbl>
## 1 exp                             26.7             0.799    46.7    6.76    1.94
## # … with abbreviated variable names ¹​`sd(msmt_temp_C, na.rm = T)`,
## #   ²​`mean(msmt_RH_percent, na.rm = T)`, ³​`sd(msmt_RH_percent, na.rm = T)`,
## #   ⁴​`mean(msmt_VPD_kPa, na.rm = T)`

Means by Day

Calculate mean values per day per tmt group.

means <- dat_no_rehab %>%
  group_by(day_n, tmt) %>%
  summarise(n_lizards = n(),
            mean_CEWL = mean(CEWL_g_m2h_mean, na.rm = T),
            sd_CEWL = sd(CEWL_g_m2h_mean, na.rm = T),
            mean_osml = mean(osmolality_mmol_kg_mean, na.rm = T),
            sd_osml = sd(osmolality_mmol_kg_mean, na.rm = T),
            mean_hct = mean(hematocrit_percent, na.rm = T),
            sd_hct = sd(hematocrit_percent, na.rm = T),
            mean_mass = mean(mass_g, na.rm = T),
            sd_mass = sd(mass_g, na.rm = T),
            mean_SMI = mean(SMI, na.rm = T),
            sd_SMI = sd(SMI, na.rm = T)) %>%
  mutate(se_CEWL = (sd_CEWL/sqrt(n_lizards)),
         se_osml = (sd_osml/sqrt(n_lizards)),
         se_hct = (sd_hct/sqrt(n_lizards)),
         se_mass = (sd_mass/sqrt(n_lizards)),
         se_SMI = (sd_SMI/sqrt(n_lizards)))
## `summarise()` has grouped output by 'day_n'. You can override using the
## `.groups` argument.
# get rid of non-defined points
# these are from when not all variables were measured for a given date
means$mean_CEWL[is.nan(means$mean_CEWL)] <- NA
means$mean_osml[is.nan(means$mean_osml)] <- NA
means$mean_hct[is.nan(means$mean_hct)] <- NA
means$mean_mass[is.nan(means$mean_mass)] <- NA
means$mean_SMI[is.nan(means$mean_SMI)] <- NA
means
## # A tibble: 24 × 18
## # Groups:   day_n [6]
##    day_n tmt      n_liz…¹ mean_…² sd_CEWL mean_…³ sd_osml mean_…⁴ sd_hct mean_…⁵
##    <dbl> <fct>      <int>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>   <dbl>
##  1     0 Cool Hu…      34    20.9    4.78    351.    20.3    39.6   5.30    11.6
##  2     0 Hot Hum…      33    21.4    4.85    347.    18.7    37.9   5.46    11.6
##  3     0 Cool Dr…      33    20.0    6.08    346.    20.6    39.3   5.96    11.8
##  4     0 Hot Dry…      34    20.9    5.93    347.    17.9    38.9   5.05    12.0
##  5     4 Cool Hu…      34    NA     NA       356.    24.5    34.5   5.34    11.1
##  6     4 Hot Hum…      33    NA     NA       359.    22.5    32.0   5.38    10.6
##  7     4 Cool Dr…      33    NA     NA       355.    27.0    35.0   7.02    11.1
##  8     4 Hot Dry…      34    NA     NA       361.    25.8    33.1   5.00    10.6
##  9     5 Cool Hu…      34    NA     NA        NA     NA      NA    NA       10.8
## 10     5 Hot Hum…      33    NA     NA        NA     NA      NA    NA       10.2
## # … with 14 more rows, 8 more variables: sd_mass <dbl>, mean_SMI <dbl>,
## #   sd_SMI <dbl>, se_CEWL <dbl>, se_osml <dbl>, se_hct <dbl>, se_mass <dbl>,
## #   se_SMI <dbl>, and abbreviated variable names ¹​n_lizards, ²​mean_CEWL,
## #   ³​mean_osml, ⁴​mean_hct, ⁵​mean_mass
# get only means for the very end
end_means <- means %>%
  dplyr::filter(day_n == 8)
write.csv(end_means, "./results_statistics/exp_end_means.csv")

End Values Only

Select all raw measurements for only day=8 values.

end_vals <- dat %>%
  dplyr::filter(day_n == 8)

delta CEWL

Get a df that only has complete observations that include CEWL values (only obs from before and after the experiment). Then, calculate the change (delta) in CEWL from before to after the experiment. Because we only measured CEWL at those two time points, it makes more sense to assess the amount of change in CEWL for each lizard, rather than measuring the change over time.

start_CEWL <- dat_no_rehab %>%
  dplyr::filter(day_n == 0) %>%
  dplyr::select(individual_ID, start_CEWL = CEWL_g_m2h_mean)
dat_no_rehab_deltaCEWL <- dat_no_rehab %>% # initiate new df
  dplyr::filter(complete.cases(CEWL_g_m2h_mean)) %>% # only use obs incl CEWL
  dplyr::filter(day_n == 8) %>% # get only obs for post-exp
  left_join(start_CEWL, by = 'individual_ID') %>% # add start CEWL to both obs for each lizard
  mutate(delta_CEWL = CEWL_g_m2h_mean - start_CEWL) # calculate deltaCEWL after-before experiment

Experiment Models

We predicted that there would be effects of day, humidity treatment, temperature treatment, and treatment VPD. However, we can’t use the standard backwards model selection because the three treatment variables are collinear (VIF much higher than acceptable) and it leads to issues with changing the sign of the estimates when all three are included together. So, we will run singular models with each treatment variable alone:

response ~ dayhumidity response ~ daytemperature response ~ day*VPD

Then, we will use AIC, RMSE, and R-sq to assess which treatment effect is most important to that response variable.

Body Mass

Building

Build each treatment effect model.

mass_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')
mass_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')
mass_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')

Assumptions

Check linear regression assumptions/conditions.

# distribution of SMI
hist(dat_no_rehab$mass_g)

# VPD model
plot(mass_mod_VPD)

simple.eda(residuals(mass_mod_VPD))

shapiro.test(residuals(mass_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mass_mod_VPD)
## W = 0.93358, p-value < 2.2e-16
# humidity model
plot(mass_mod_hum)

simple.eda(residuals(mass_mod_hum))

shapiro.test(residuals(mass_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mass_mod_hum)
## W = 0.93039, p-value < 2.2e-16
# temperature model
plot(mass_mod_temp)

simple.eda(residuals(mass_mod_temp))

shapiro.test(residuals(mass_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mass_mod_temp)
## W = 0.90802, p-value < 2.2e-16

Normality is violated, but linearity, equal error variance, and independence are all good.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

We…

  • calculate RMSE manually
  • use the r.squaredGLMM function in the MuMIn package to get the marginal R^2, which is how much of the total variance is explained by fixed effects
  • use the aictab function in the AICmodavg package to get AIC and deltaAIC values
  • get the sum of squares, F, and p-values for each variable from the anova table for each model
# calculate RMSE & R^2
mass_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(mass_mod_VPD))^2)),
                                    sqrt(mean((residuals(mass_mod_hum))^2)),
                                    sqrt(mean((residuals(mass_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(mass_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(mass_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(mass_mod_temp)[,"R2m"]))
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
# calculate AIC
mass_models <- list(mass_mod_VPD, mass_mod_hum, mass_mod_temp)
EXP_mod_names <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ))
mass_AICc <- data.frame(aictab(cand.set = mass_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = mass_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
mass_across <- mass_RMSE_Rsq %>%
  left_join(mass_AICc, by = 'model') %>%
  mutate(response = "Body Mass (g)") %>%
  arrange(Delta_AICc)

# calculate F & p-values
mass_VPD_p <- data.frame(anova(mass_mod_VPD, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
mass_hum_p <- data.frame(anova(mass_mod_hum, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
mass_temp_p <- data.frame(anova(mass_mod_temp, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
mass_within <- mass_VPD_p %>%
  rbind(mass_hum_p) %>%
  rbind(mass_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Body Mass (g)")

Body Condition

Building

Build each treatment effect model.

SMI_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
SMI_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
SMI_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))

Assumptions

Check linear regression assumptions/conditions.

# distribution of SMI
hist(dat_no_rehab$SMI)

# VPD model
plot(SMI_mod_VPD)

simple.eda(residuals(SMI_mod_VPD))

shapiro.test(residuals(SMI_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(SMI_mod_VPD)
## W = 0.96052, p-value = 6.691e-14
# humidity model
plot(SMI_mod_hum)

simple.eda(residuals(SMI_mod_hum))

shapiro.test(residuals(SMI_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(SMI_mod_hum)
## W = 0.95569, p-value = 7.603e-15
# temperature model
plot(SMI_mod_temp)

simple.eda(residuals(SMI_mod_temp))

shapiro.test(residuals(SMI_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(SMI_mod_temp)
## W = 0.93877, p-value < 2.2e-16

Normality is violated, but linearity, equal error variance, and independence are all good.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

We…

  • calculate RMSE manually
  • use the r.squaredGLMM function in the MuMIn package to get the marginal R^2, which is how much of the total variance is explained by fixed effects
  • use the aictab function in the AICmodavg package to get AIC and deltaAIC values
  • get the sum of squares, F, and p-values for each variable from the anova table for each model
# calculate RMSE & R^2
SMI_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(SMI_mod_VPD))^2)),
                                    sqrt(mean((residuals(SMI_mod_hum))^2)),
                                    sqrt(mean((residuals(SMI_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(SMI_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(SMI_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(SMI_mod_temp)[,"R2m"]))

# calculate AIC
SMI_models <- list(SMI_mod_VPD, SMI_mod_hum, SMI_mod_temp)
EXP_mod_names <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ))
SMI_AICc <- data.frame(aictab(cand.set = SMI_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = SMI_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
SMI_across <- SMI_RMSE_Rsq %>%
  left_join(SMI_AICc, by = 'model') %>%
  mutate(response = "Body Condition (g')") %>%
  arrange(Delta_AICc)

# calculate F & p-values
SMI_VPD_p <- data.frame(anova(SMI_mod_VPD, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
SMI_hum_p <- data.frame(anova(SMI_mod_hum, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
SMI_temp_p <- data.frame(anova(SMI_mod_temp, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
SMI_within <- SMI_VPD_p %>%
  rbind(SMI_hum_p) %>%
  rbind(SMI_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Body Condition (g')")

Hematocrit

Building

Build each treatment effect model.

hct_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                               hematocrit_percent ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
hct_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                              hematocrit_percent ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
hct_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                              hematocrit_percent ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))

Assumptions

Check linear regression assumptions/conditions.

# distribution of 
hist(dat_no_rehab$hematocrit_percent)

# VPD model
plot(hct_mod_VPD)

simple.eda(residuals(hct_mod_VPD))

shapiro.test(residuals(hct_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_mod_VPD)
## W = 0.99651, p-value = 0.542
# humidity model
plot(hct_mod_hum)

simple.eda(residuals(hct_mod_hum))

shapiro.test(residuals(hct_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_mod_hum)
## W = 0.99619, p-value = 0.4584
# temperature model
plot(hct_mod_temp)

simple.eda(residuals(hct_mod_temp))

shapiro.test(residuals(hct_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_mod_temp)
## W = 0.99478, p-value = 0.1975

All assumptions, normality, linearity, equal error variance, and independence are all good.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

# calculate RMSE & R^2
hct_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(hct_mod_VPD))^2)),
                                    sqrt(mean((residuals(hct_mod_hum))^2)),
                                    sqrt(mean((residuals(hct_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(hct_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(hct_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(hct_mod_temp)[,"R2m"]))

# calculate AIC
hct_models <- list(hct_mod_VPD, hct_mod_hum, hct_mod_temp)
hct_AICc <- data.frame(aictab(cand.set = hct_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = hct_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
hct_across <- hct_RMSE_Rsq %>%
  left_join(hct_AICc, by = 'model') %>%
  mutate(response = "Hematocrit (%)") %>%
  arrange(Delta_AICc)

# calculate F & p-values
hct_VPD_p <- data.frame(anova(hct_mod_VPD, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
hct_hum_p <- data.frame(anova(hct_mod_hum, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
hct_temp_p <- data.frame(anova(hct_mod_temp, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
hct_within <- hct_VPD_p %>%
  rbind(hct_hum_p) %>%
  rbind(hct_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Hematocrit (%)")

Osmolality

Building

Build each treatment effect model.

osml_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                               osmolality_mmol_kg_mean ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
osml_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                               osmolality_mmol_kg_mean ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
osml_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                               osmolality_mmol_kg_mean ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))

Assumptions

Check linear regression assumptions/conditions.

# distribution of 
hist(dat_no_rehab$osmolality_mmol_kg_mean)

# VPD model
plot(osml_mod_VPD)

simple.eda(residuals(osml_mod_VPD))

shapiro.test(residuals(osml_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(osml_mod_VPD)
## W = 0.9769, p-value = 6.836e-06
# humidity model
plot(osml_mod_hum)

simple.eda(residuals(osml_mod_hum))

shapiro.test(residuals(osml_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(osml_mod_hum)
## W = 0.97746, p-value = 8.914e-06
# temperature model
plot(osml_mod_temp)

simple.eda(residuals(osml_mod_temp))

shapiro.test(residuals(osml_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(osml_mod_temp)
## W = 0.97346, p-value = 1.44e-06

Normality is violated, but linearity, equal error variance, and independence are all okay.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

# calculate RMSE & R^2
osml_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(osml_mod_VPD))^2)),
                                    sqrt(mean((residuals(osml_mod_hum))^2)),
                                    sqrt(mean((residuals(osml_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(osml_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(osml_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(osml_mod_temp)[,"R2m"]))

# calculate AIC
osml_models <- list(osml_mod_VPD, osml_mod_hum, osml_mod_temp)
osml_AICc <- data.frame(aictab(cand.set = osml_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = osml_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
osml_across <- osml_RMSE_Rsq %>%
  left_join(osml_AICc, by = 'model') %>%
  mutate(response = "Plasma Osmolality (mmol/kg)") %>%
  arrange(Delta_AICc)

# calculate F & p-values
osml_VPD_p <- data.frame(anova(osml_mod_VPD, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
osml_hum_p <- data.frame(anova(osml_mod_hum, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
osml_temp_p <- data.frame(anova(osml_mod_temp, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
osml_within <- osml_VPD_p %>%
  rbind(osml_hum_p) %>%
  rbind(osml_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Plasma Osmolality (mmol/kg)")

Body Temp

I need to double check whether CEWL has a relationship with body temperature at the point of measurement.

ggplot(dat_no_rehab) + 
  geom_point(aes(x = CEWL_g_m2h_mean,
                 y = cloacal_temp_C,
                 color = day_n))
## Warning: Removed 537 rows containing missing values (`geom_point()`).

Test an lm of raw CEWL ~ body temp for day 8 measurements only.

body_temp_test_dat <- end_vals %>%
  dplyr::filter(complete.cases(CEWL_g_m2h_mean))

CEWL_body_temp_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                         CEWL_g_m2h_mean ~ cloacal_temp_C * tmt +
                           (1|trial_number))
summary(CEWL_body_temp_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ cloacal_temp_C * tmt + (1 | trial_number)
##    Data: body_temp_test_dat
## 
## REML criterion at convergence: 839.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.02426 -0.58867 -0.09161  0.46172  3.09926 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  4.977   2.231   
##  Residual                 36.926   6.077   
## Number of obs: 133, groups:  trial_number, 5
## 
## Fixed effects:
##                                        Estimate Std. Error        df t value
## (Intercept)                            42.83759   26.65061 123.21107   1.607
## cloacal_temp_C                         -0.49420    1.05450 122.93003  -0.469
## tmtHot Humid (1.1 kPa)                -16.46275   34.17703 122.09739  -0.482
## tmtCool Dry (2.5 kPa)                 -23.23095   47.80497 123.10602  -0.486
## tmtHot Dry (3.8 kPa)                   -3.95989   36.73317 123.70420  -0.108
## cloacal_temp_C:tmtHot Humid (1.1 kPa)   0.90966    1.35821 122.13178   0.670
## cloacal_temp_C:tmtCool Dry (2.5 kPa)    0.65396    1.86629 123.06213   0.350
## cloacal_temp_C:tmtHot Dry (3.8 kPa)    -0.06271    1.44351 123.65911  -0.043
##                                       Pr(>|t|)
## (Intercept)                              0.111
## cloacal_temp_C                           0.640
## tmtHot Humid (1.1 kPa)                   0.631
## tmtCool Dry (2.5 kPa)                    0.628
## tmtHot Dry (3.8 kPa)                     0.914
## cloacal_temp_C:tmtHot Humid (1.1 kPa)    0.504
## cloacal_temp_C:tmtCool Dry (2.5 kPa)     0.727
## cloacal_temp_C:tmtHot Dry (3.8 kPa)      0.965
## 
## Correlation of Fixed Effects:
##             (Intr) clc__C tHH(1k tCD(2k tHD(3k c__H(k c__C:CD(k
## clocl_tmp_C -0.999                                             
## tmHH(1.1kP) -0.745  0.744                                      
## tmCD(2.5kP) -0.528  0.528  0.434                               
## tmHD(3.8kP) -0.694  0.694  0.564  0.422                        
## c__C:HH(1.k  0.741 -0.741 -0.999 -0.432 -0.562                 
## c__C:CD(2.k  0.536 -0.536 -0.439 -0.999 -0.427  0.438          
## c__C:HD(3.k  0.699 -0.700 -0.567 -0.425 -0.999  0.566  0.430
anova(CEWL_body_temp_mod, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                     Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)    
## cloacal_temp_C      230.75  230.75     1 115.89  6.2491  0.01383 *  
## tmt                3115.05 1038.35     3 122.21 28.1193 6.75e-14 ***
## cloacal_temp_C:tmt   27.57    9.19     3 122.56  0.2489  0.86199    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(body_temp_test_dat) + 
  geom_point(aes(x = msmt_temp_C,
                 y = cloacal_temp_C)) +
  xlim(23, 28) + ylim(23,28)

ggplot(body_temp_test_dat) + 
  geom_point(aes(x = msmt_temp_C,
                 y = CEWL_g_m2h_mean,
                 color = tmt))

Do the same exact stats for ambient temp and VPD at the time of measurement:

CEWL_msmt_temp_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                                     CEWL_g_m2h_mean ~ msmt_temp_C * tmt +
                                       (1|trial_number))
CEWL_msmt_VPD_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                                     CEWL_g_m2h_mean ~ msmt_VPD_kPa * tmt +
                                       (1|trial_number))

summary(CEWL_msmt_temp_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ msmt_temp_C * tmt + (1 | trial_number)
##    Data: body_temp_test_dat
## 
## REML criterion at convergence: 824.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8115 -0.5673 -0.0847  0.4720  3.5376 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  7.298   2.701   
##  Residual                 33.499   5.788   
## Number of obs: 133, groups:  trial_number, 5
## 
## Fixed effects:
##                                    Estimate Std. Error       df t value
## (Intercept)                        -39.6305    48.2297 112.4351  -0.822
## msmt_temp_C                          2.6733     1.8380 112.5242   1.454
## tmtHot Humid (1.1 kPa)             -73.4611    65.2521 121.4757  -1.126
## tmtCool Dry (2.5 kPa)              -18.6586    60.6146 121.3298  -0.308
## tmtHot Dry (3.8 kPa)               -61.0925    64.0861 121.2408  -0.953
## msmt_temp_C:tmtHot Humid (1.1 kPa)   3.0314     2.4837 121.4794   1.221
## msmt_temp_C:tmtCool Dry (2.5 kPa)    0.4532     2.3073 121.3153   0.196
## msmt_temp_C:tmtHot Dry (3.8 kPa)     2.1152     2.4428 121.2415   0.866
##                                    Pr(>|t|)
## (Intercept)                           0.413
## msmt_temp_C                           0.149
## tmtHot Humid (1.1 kPa)                0.262
## tmtCool Dry (2.5 kPa)                 0.759
## tmtHot Dry (3.8 kPa)                  0.342
## msmt_temp_C:tmtHot Humid (1.1 kPa)    0.225
## msmt_temp_C:tmtCool Dry (2.5 kPa)     0.845
## msmt_temp_C:tmtHot Dry (3.8 kPa)      0.388
## 
## Correlation of Fixed Effects:
##             (Intr) msm__C tHH(1k tCD(2k tHD(3k m__H(k m__C:CD(k
## msmt_temp_C -0.999                                             
## tmHH(1.1kP) -0.593  0.592                                      
## tmCD(2.5kP) -0.614  0.613  0.467                               
## tmHD(3.8kP) -0.580  0.580  0.445  0.478                        
## m__C:HH(1.k  0.594 -0.594 -1.000 -0.467 -0.445                 
## m__C:CD(2.k  0.614 -0.614 -0.467 -1.000 -0.478  0.467          
## m__C:HD(3.k  0.579 -0.579 -0.444 -0.478 -1.000  0.445  0.478
summary(CEWL_msmt_VPD_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ msmt_VPD_kPa * tmt + (1 | trial_number)
##    Data: body_temp_test_dat
## 
## REML criterion at convergence: 817.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8187 -0.5552 -0.0438  0.3596  3.2064 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 11.75    3.428   
##  Residual                 34.18    5.846   
## Number of obs: 133, groups:  trial_number, 5
## 
## Fixed effects:
##                                     Estimate Std. Error      df t value
## (Intercept)                           21.583     14.417  63.780   1.497
## msmt_VPD_kPa                           4.952      8.054  67.896   0.615
## tmtHot Humid (1.1 kPa)               -24.082     16.730 120.618  -1.439
## tmtCool Dry (2.5 kPa)                 -9.284     16.225 120.721  -0.572
## tmtHot Dry (3.8 kPa)                 -25.308     16.211 120.503  -1.561
## msmt_VPD_kPa:tmtHot Humid (1.1 kPa)   17.094      9.373 120.640   1.824
## msmt_VPD_kPa:tmtCool Dry (2.5 kPa)     1.504      9.025 120.677   0.167
## msmt_VPD_kPa:tmtHot Dry (3.8 kPa)     11.097      9.093 120.529   1.220
##                                     Pr(>|t|)  
## (Intercept)                           0.1393  
## msmt_VPD_kPa                          0.5407  
## tmtHot Humid (1.1 kPa)                0.1526  
## tmtCool Dry (2.5 kPa)                 0.5682  
## tmtHot Dry (3.8 kPa)                  0.1211  
## msmt_VPD_kPa:tmtHot Humid (1.1 kPa)   0.0707 .
## msmt_VPD_kPa:tmtCool Dry (2.5 kPa)    0.8679  
## msmt_VPD_kPa:tmtHot Dry (3.8 kPa)     0.2247  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) m_VPD_ tHH(1k tCD(2k tHD(3k m_VH(k m_VPD_P:CD(k
## msmt_VPD_kP  -0.992                                                
## tmHH(1.1kP)  -0.538  0.538                                         
## tmCD(2.5kP)  -0.524  0.523  0.472                                  
## tmHD(3.8kP)  -0.535  0.535  0.476  0.490                           
## m_VPD_P:H(k   0.538 -0.542 -0.996 -0.470 -0.474                    
## m_VPD_P:CD(k  0.527 -0.531 -0.474 -0.996 -0.492  0.476             
## m_VPD_P:HD(k  0.529 -0.533 -0.473 -0.488 -0.996  0.475  0.493
anova(CEWL_msmt_temp_mod, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                 Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## msmt_temp_C      320.3  320.27     1  59.349  9.5605  0.003027 ** 
## tmt             3464.5 1154.85     3 121.436 34.4742 3.438e-16 ***
## msmt_temp_C:tmt   65.4   21.78     3 121.538  0.6503  0.584244    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(CEWL_msmt_VPD_mod, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                  Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## msmt_VPD_kPa       91.2   91.22     1  32.035  2.6689    0.1121    
## tmt              3536.4 1178.81     3 121.547 34.4879 3.372e-16 ***
## msmt_VPD_kPa:tmt  152.9   50.97     3 121.397  1.4911    0.2204    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CEWL

First, get a separate df of end-experiment delta CEWL for each tmt group:

dat_no_rehab_deltaCEWL_post <- dat_no_rehab_deltaCEWL %>%
  dplyr::filter(day_n == 8)
HH8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Hu")
HD8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Dr")
CH8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool H")
CD8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool D")

t-tests

I use t-tests to assess whether the change in CEWL from the experiment is significantly different from zero for each treatment group.

CEWL_ttest_HH <- t.test(HH8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_HH
## 
##  One Sample t-test
## 
## data:  HH8$delta_CEWL
## t = 8.7341, df = 32, p-value = 5.573e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  11.76821 18.92672
## sample estimates:
## mean of x 
##  15.34746
CEWL_ttest_HD <- t.test(HD8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_HD
## 
##  One Sample t-test
## 
## data:  HD8$delta_CEWL
## t = 2.513, df = 33, p-value = 0.01703
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.7042497 6.6929758
## sample estimates:
## mean of x 
##  3.698613
CEWL_ttest_CH <- t.test(CH8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_CH
## 
##  One Sample t-test
## 
## data:  CH8$delta_CEWL
## t = 8.7488, df = 32, p-value = 5.363e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   7.296566 11.725292
## sample estimates:
## mean of x 
##  9.510929
CEWL_ttest_CD <- t.test(CD8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_CD
## 
##  One Sample t-test
## 
## data:  CD8$delta_CEWL
## t = 2.75, df = 32, p-value = 0.009721
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.9302207 6.2446177
## sample estimates:
## mean of x 
##  3.587419

Building

Build each treatment effect model. These are equivalent to the models built for body condition, hematocrit, and plasma osmolality, but the effect of time is not included because we assess CEWL as delta CEWL.

CEWL_mod_VPD <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ VPD_kPa +
                              (1|trial_number))
CEWL_mod_hum <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ humidity_tmt +
                              (1|trial_number))
CEWL_mod_temp <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ temp_tmt +
                              (1|trial_number))

Assumptions

Check linear regression assumptions/conditions.

# distribution of 
hist(dat_no_rehab_deltaCEWL$delta_CEWL)

# VPD model
plot(CEWL_mod_VPD)

simple.eda(residuals(CEWL_mod_VPD))

shapiro.test(residuals(CEWL_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_mod_VPD)
## W = 0.97332, p-value = 0.01016
# humidity model
plot(CEWL_mod_hum)

simple.eda(residuals(CEWL_mod_hum))

shapiro.test(residuals(CEWL_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_mod_hum)
## W = 0.98319, p-value = 0.1002
# temperature model
plot(CEWL_mod_temp)

simple.eda(residuals(CEWL_mod_temp))

shapiro.test(residuals(CEWL_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_mod_temp)
## W = 0.98262, p-value = 0.08761

All assumptions are fine.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

# calculate RMSE & R^2
CEWL_RMSE_Rsq <- data.frame(model = 
                               c('VPD',
                                 'Humidity',
                                 'Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(CEWL_mod_VPD))^2)),
                                    sqrt(mean((residuals(CEWL_mod_hum))^2)),
                                    sqrt(mean((residuals(CEWL_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(CEWL_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(CEWL_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(CEWL_mod_temp)[,"R2m"]))

# calculate AIC
CEWL_models <- list(CEWL_mod_VPD, CEWL_mod_hum, CEWL_mod_temp)
CEWL_mod_names <- data.frame(model = 
                               c('VPD',
                                 'Humidity',
                                 'Temp'
                                 ))
CEWL_AICc <- data.frame(aictab(cand.set = CEWL_models, 
                                 modnames = CEWL_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = CEWL_models, modnames = CEWL_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
CEWL_across <- CEWL_RMSE_Rsq %>%
  left_join(CEWL_AICc, by = 'model') %>%
  mutate(response = "deltaCEWL") %>%
  arrange(Delta_AICc)

# calculate F & p-values
CEWL_VPD_p <- data.frame(anova(CEWL_mod_VPD, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'VPD',
         predictor = rownames(.))
CEWL_hum_p <- data.frame(anova(CEWL_mod_hum, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Humidity',
         predictor = rownames(.))
CEWL_temp_p <- data.frame(anova(CEWL_mod_temp, 
                              type = "1", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Temp',
         predictor = rownames(.))

# save within model values
CEWL_within <- CEWL_VPD_p %>%
  rbind(CEWL_hum_p) %>%
  rbind(CEWL_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "deltaCEWL")

Group Export

Put all the model statistics into one df/csv - one for among-model comparisons, and one for within-model stats.

# comparisons of all models for each variable
experiment_model_compare <- CEWL_across %>%
  rbind(osml_across) %>%
  rbind(hct_across) %>%
  rbind(mass_across) %>%
  rbind(SMI_across) %>%
  dplyr::select(response, model, 
                RMSE, Rsq, AICc, Delta_AICc) %>%
  mutate(RMSE = round(RMSE, 2), 
         Rsq = round(Rsq, 2), 
         AICc = round(AICc, 2), 
         Delta_AICc = round(Delta_AICc, 2)) 
write.csv(experiment_model_compare, 
          "./results_statistics/exp_model_comparisons.csv")

# all of the actual model results
experiment_model_values <- CEWL_within %>%
  rbind(osml_within) %>%
  rbind(hct_within) %>%
  rbind(mass_within) %>%
  rbind(SMI_within) %>%
  dplyr::select(response, model, predictor, 
                seq_sum_of_squares = Sum.Sq, 
                df, 
                F_statistic = F.value, 
                p_value = Pr..F.) %>%
  mutate(seq_sum_of_squares = round(seq_sum_of_squares, 0), 
         F_statistic = round(F_statistic, 2)) 
write.csv(experiment_model_values, "./results_statistics/exp_model_values.csv")

Effect Estimates

End Value CIs

Now, we can use the emmeans function from the emmeans package to estimate marginal means and confidence intervals for the values among treatment groups at the end of the experiment. But, to get for each treatment group, we need to run a new model with day 8 data only and tmt as a single, 4-category variable.

# Body Mass
mass_mod_end <- lmerTest::lmer(data = end_vals,
                              mass_g ~ tmt + 
                              (1|trial_number))
## boundary (singular) fit: see help('isSingular')
mass_emmeans <- data.frame(emmeans(mass_mod_end, "tmt")) %>%
  mutate(response = "Body Mass (g)")
mass_pairwise <- data.frame(pairs(emmeans(mass_mod_end, "tmt"))) %>%
  mutate(response = "Body Mass (g)")

# Body Condition
SMI_mod_end <- lmerTest::lmer(data = end_vals,
                              SMI ~ tmt + 
                              (1|trial_number))
SMI_emmeans <- data.frame(emmeans(SMI_mod_end, "tmt")) %>%
  mutate(response = "Body Condition (g')")
SMI_pairwise <- data.frame(pairs(emmeans(SMI_mod_end, "tmt"))) %>%
  mutate(response = "Body Condition (g')")

# Hematocrit
hct_mod_end <- lmerTest::lmer(data = end_vals,
                              hematocrit_percent ~ tmt + 
                              (1|trial_number))
hct_emmeans <- data.frame(emmeans(hct_mod_end, "tmt")) %>%
  mutate(response = "Hematocrit (%)")
hct_pairwise <- data.frame(pairs(emmeans(hct_mod_end, "tmt"))) %>%
  mutate(response = "Hematocrit (%)")

# Plasma Osmolality
osml_mod_end <- lmerTest::lmer(data = end_vals,
                              osmolality_mmol_kg_mean ~ tmt + 
                              (1|trial_number))
osml_emmeans <- data.frame(emmeans(osml_mod_end, "tmt")) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")
osml_pairwise <- data.frame(pairs(emmeans(osml_mod_end, "tmt"))) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")

# CEWL
CEWL_mod_end <- lmerTest::lmer(data = end_vals,
                              CEWL_g_m2h_mean ~ tmt + 
                              (1|trial_number))
CEWL_emmeans <- data.frame(emmeans(CEWL_mod_end, "tmt")) %>%
  mutate(response = "CEWL (g/m2h)")
CEWL_pairwise <- data.frame(pairs(emmeans(CEWL_mod_end, "tmt"))) %>%
  mutate(response = "CEWL (g/m2h)")

Rate of Change

# Body Condition
SMI_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
SMI_emtrends <- data.frame(emtrends(SMI_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Body Condition (g')")
SMI_pairtrend <- data.frame(pairs(emtrends(SMI_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Body Condition (g')")

# Hematocrit
hct_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              hematocrit_percent ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
hct_emtrends <- data.frame(emtrends(hct_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Hematocrit (%)")
hct_pairtrend <- data.frame(pairs(emtrends(hct_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Hematocrit (%)")

# Plasma Osmolality
osml_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              osmolality_mmol_kg_mean ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
osml_emtrends <- data.frame(emtrends(osml_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")
osml_pairtrend <- data.frame(pairs(emtrends(osml_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")

# CEWL
CEWL_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              CEWL_g_m2h_mean ~ day_n * tmt + 
                              (1|trial_number))
CEWL_emtrends <- data.frame(emtrends(CEWL_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "CEWL (g/m2h)")
CEWL_pairtrend <- data.frame(pairs(emtrends(CEWL_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "CEWL (g/m2h)")

# put together
all_emtrends <- CEWL_emtrends %>%
  rbind(osml_emtrends) %>%
  rbind(hct_emtrends) %>%
  rbind(SMI_emtrends) %>%
  mutate(confint95 = paste(round(lower.CL, 2), round(upper.CL, 2), sep = ", ")) %>%
  dplyr::select(response, tmt, 
                per_day_change = day_n.trend,
                confint95,
                SE, df)
write.csv(all_emtrends, "./results_statistics/exp_emtrends_per_day_change.csv")

Given the daily trends for plasma osmolality, in total for acclimation change in osml had 1.04 (0.13x8) to 13 (1.587x8) mmol/kg of change.

CEWL ~ osmolality

Let’s compare the relationship of CEWL ~ osmolality post-experiment to what we found pre-experiment in the “analysis_capture” file. The treatment-specific sub-dataframes were created in the “Body Temp” model subsection above.

CEWL_osml <- lmerTest::lmer(data = end_vals, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: end_vals
## 
## REML criterion at convergence: 855.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1565 -0.6405 -0.2787  0.6163  3.5785 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  5.66    2.379   
##  Residual                 58.25    7.632   
## Number of obs: 123, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)             18.53375   10.95312 66.51336   1.692   0.0953 .
## osmolality_mmol_kg_mean  0.02904    0.03081 72.84361   0.942   0.3491  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.993
anova(CEWL_osml, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 45.195  45.195     1 69.109  0.7759 0.3814
HH_CEWL_osml <- lmerTest::lmer(data = HH8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
## boundary (singular) fit: see help('isSingular')
summary(HH_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: HH8
## 
## REML criterion at convergence: 222.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1357 -0.6398 -0.1682  0.4261  2.4440 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  0.00    0.00    
##  Residual                 62.57    7.91    
## Number of obs: 32, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)             19.77938   23.68747 30.00000   0.835    0.410
## osmolality_mmol_kg_mean  0.04658    0.06734 30.00000   0.692    0.494
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.998
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(HH_CEWL_osml, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 24.554  24.554     1 18.124  0.3924 0.5388
HD_CEWL_osml <- lmerTest::lmer(data = HD8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(HD_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: HD8
## 
## REML criterion at convergence: 171.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55684 -0.59339 -0.07624  0.40035  2.96626 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  6.99    2.644   
##  Residual                 11.28    3.358   
## Number of obs: 31, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)              2.11489   11.95116 25.69290   0.177    0.861  
## osmolality_mmol_kg_mean  0.06101    0.03291 26.24328   1.854    0.075 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.994
anova(HD_CEWL_osml, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 31.917  31.917     1 26.042  2.8298 0.1045
CH_CEWL_osml <- lmerTest::lmer(data = CH8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(CH_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: CH8
## 
## REML criterion at convergence: 192.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2869 -0.7878 -0.1247  0.6553  2.0089 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  6.866   2.620   
##  Residual                 30.995   5.567   
## Number of obs: 30, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)             22.02423   13.70141 27.96594   1.607    0.119
## osmolality_mmol_kg_mean  0.02489    0.03901 27.97274   0.638    0.529
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.993
anova(CH_CEWL_osml, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 11.376  11.376     1 27.97   0.367 0.5495
CD_CEWL_osml <- lmerTest::lmer(data = CD8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(CD_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: CD8
## 
## REML criterion at convergence: 170
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8562 -0.5019  0.1381  0.4084  1.8415 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 21.14    4.598   
##  Residual                 11.34    3.368   
## Number of obs: 30, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)             -1.45589   11.16955 27.92145  -0.130   0.8972  
## osmolality_mmol_kg_mean  0.07357    0.03133 27.93722   2.348   0.0262 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.981
anova(CD_CEWL_osml, type = "1", ddf = "Kenward-Roger")
## Type I Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## osmolality_mmol_kg_mean 55.393  55.393     1 27.937  4.8837 0.03548 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Even though only data for one treatment group was significant, all have a positive relationship

CEWL ~ VPD

We also want to assess whether the vapor pressure deficit lizards experienced during acclimation has an effect on the change in CEWL.

CEWL_VPD_lm <- lm(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ VPD_kPa_tmttrial)
plot(CEWL_VPD_lm)

simple.eda(residuals(CEWL_VPD_lm))

shapiro.test(residuals(CEWL_VPD_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_VPD_lm)
## W = 0.96694, p-value = 0.002527
summary(CEWL_VPD_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ VPD_kPa_tmttrial, data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8069  -5.6289  -0.7159   5.5303  27.3517 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        13.960      1.415   9.866  < 2e-16 ***
## VPD_kPa_tmttrial   -2.959      0.594  -4.982 1.95e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.726 on 131 degrees of freedom
## Multiple R-squared:  0.1593, Adjusted R-squared:  0.1529 
## F-statistic: 24.82 on 1 and 131 DF,  p-value: 1.952e-06

Even though the data is slightly nonlinear, a linear model does a fine job explaining the data.

To be extra sure, also check for high leverage and influential points using H-values and Cook’s distance:

# compute values for observations 
high_leverage <- data.frame(H = hatvalues(CEWL_VPD_lm)) %>% 
  mutate(row = row_number())

# compute cutoff value 
h_bar <- (3*sum(high_leverage$H))/nrow(high_leverage)

# add to original dataframe 
# see which observations have extremely high leverage (if any)
high_leverage_dat <- dat_no_rehab_deltaCEWL %>%
  mutate(row = row_number()) %>%
  left_join(., high_leverage, by = "row") %>%
  dplyr::filter(H > h_bar) 
high_leverage_dat
## # A tibble: 0 × 33
## # Groups:   individual_ID [0]
## # … with 33 variables: measurement_date <date>, type <fct>,
## #   individual_ID <fct>, mass_g <dbl>, hematocrit_percent <int>,
## #   trial_number <fct>, temp_tmt <fct>, humidity_tmt <fct>, SVL_mm <int>,
## #   tmt <fct>, day_n <dbl>, day_factor <fct>, osmolality_mmol_kg_mean <dbl>,
## #   CEWL_g_m2h_mean <dbl>, msmt_temp_C <dbl>, msmt_RH_percent <dbl>,
## #   cloacal_temp_C <dbl>, msmt_temp_K <dbl>, e_s_kPa_m <dbl>, e_a_kPa_m <dbl>,
## #   msmt_VPD_kPa <dbl>, SMI <dbl>, temp_mean_tmttrial <dbl>, …
# get Cook's distance 
cooks <- data.frame(c = cooks.distance(CEWL_VPD_lm)) %>% 
  mutate(row = row_number())

# add to original dataframe 
influential <- dat_no_rehab_deltaCEWL %>%
  mutate(row = row_number()) %>% 
  left_join(., cooks, by = "row")

# see moderately influential points 
cook_mod_inf <- influential %>% 
  dplyr::filter(c>0.5) 
cook_mod_inf
## # A tibble: 0 × 33
## # Groups:   individual_ID [0]
## # … with 33 variables: measurement_date <date>, type <fct>,
## #   individual_ID <fct>, mass_g <dbl>, hematocrit_percent <int>,
## #   trial_number <fct>, temp_tmt <fct>, humidity_tmt <fct>, SVL_mm <int>,
## #   tmt <fct>, day_n <dbl>, day_factor <fct>, osmolality_mmol_kg_mean <dbl>,
## #   CEWL_g_m2h_mean <dbl>, msmt_temp_C <dbl>, msmt_RH_percent <dbl>,
## #   cloacal_temp_C <dbl>, msmt_temp_K <dbl>, e_s_kPa_m <dbl>, e_a_kPa_m <dbl>,
## #   msmt_VPD_kPa <dbl>, SMI <dbl>, temp_mean_tmttrial <dbl>, …

NO high leverage or influential points! :)

We will double check a comparison of a polynomial model, just to be sure:

CEWL_VPD_poly <- lm(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ poly(VPD_kPa_tmttrial, 2))
plot(CEWL_VPD_poly)

simple.eda(residuals(CEWL_VPD_poly))

shapiro.test(residuals(CEWL_VPD_poly))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_VPD_poly)
## W = 0.97276, p-value = 0.00896
summary(CEWL_VPD_poly)
## 
## Call:
## lm(formula = delta_CEWL ~ poly(VPD_kPa_tmttrial, 2), data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.222  -5.281  -0.914   4.981  26.576 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   8.004      0.752  10.643  < 2e-16 ***
## poly(VPD_kPa_tmttrial, 2)1  -43.470      8.672  -5.013 1.72e-06 ***
## poly(VPD_kPa_tmttrial, 2)2  -14.069      8.672  -1.622    0.107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.672 on 130 degrees of freedom
## Multiple R-squared:  0.176,  Adjusted R-squared:  0.1633 
## F-statistic: 13.88 on 2 and 130 DF,  p-value: 3.442e-06

LINE assumptions are equally-well-met.

The polynomial factor is not significant, but the R-sq is slightly higher for the poly model.

Compare RMSE and AIC:

sqrt(mean((residuals(CEWL_VPD_lm))^2))
## [1] 8.660173
sqrt(mean((residuals(CEWL_VPD_poly))^2))
## [1] 8.573822
CEWL_VPD_models <- list(CEWL_VPD_lm, CEWL_VPD_poly)
CEWL_VPD_mod_names <- data.frame(model = 
                               c('linear',
                                 'polynomial'
                                 ))
CEWL_VPD_AICc <- data.frame(aictab(cand.set = CEWL_VPD_models, 
                                 modnames = CEWL_VPD_mod_names))
CEWL_VPD_AICc
##        model K     AICc Delta_AICc  ModelLik    AICcWt        LL    Cum.Wt
## 2 polynomial 4 957.3080  0.0000000 1.0000000 0.5669881 -474.4977 0.5669881
## 1     linear 3 957.8471  0.5391459 0.7637056 0.4330119 -475.8305 1.0000000

RMSE is slightly lower for the polynomial model. But, AIC is not meaningfully different between the two versions.

I’ll use the lm as the best model/

summary(CEWL_VPD_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ VPD_kPa_tmttrial, data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8069  -5.6289  -0.7159   5.5303  27.3517 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        13.960      1.415   9.866  < 2e-16 ***
## VPD_kPa_tmttrial   -2.959      0.594  -4.982 1.95e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.726 on 131 degrees of freedom
## Multiple R-squared:  0.1593, Adjusted R-squared:  0.1529 
## F-statistic: 24.82 on 1 and 131 DF,  p-value: 1.952e-06
anova(CEWL_VPD_lm)
## Analysis of Variance Table
## 
## Response: delta_CEWL
##                   Df Sum Sq Mean Sq F value    Pr(>F)    
## VPD_kPa_tmttrial   1 1889.7 1889.67  24.817 1.952e-06 ***
## Residuals        131 9974.8   76.14                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(estimate = -2.96, SE = 0.6, df = (1, 131), F = 24.82, p < 0.0001, Rsq = 0.16)

Recovery Models

I want to know how the 2-day recovery period affects physiology relative to post- and pre- experiment. To do this, I’ll use two-sided t-tests comparing recovery delta values to the hypothesis of mu=0.

SMI

SMI_rmod_post_exp <- t.test(recovery_v_post_exp$delta_SMI_10_8, 
                           mu = 0, alternative = "two.sided")
SMI_rmod_post_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_post_exp$delta_SMI_10_8
## t = 6.677, df = 131, p-value = 6.292e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3171205 0.5841444
## sample estimates:
## mean of x 
## 0.4506324
SMI_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_SMI_10_0, 
                           mu = 0, alternative = "two.sided")
SMI_rmod_pre_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_pre_exp$delta_SMI_10_0
## t = -11.542, df = 131, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.3426322 -0.9497445
## sample estimates:
## mean of x 
## -1.146188

Hematocrit

hct_rmod_post_exp <- t.test(recovery_v_post_exp$delta_hct_10_8, 
                           mu = 0, alternative = "two.sided")
hct_rmod_post_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_post_exp$delta_hct_10_8
## t = -4.2083, df = 126, p-value = 4.85e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -3.033119 -1.092866
## sample estimates:
## mean of x 
## -2.062992
hct_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_hct_10_0, 
                           mu = 0, alternative = "two.sided")
hct_rmod_pre_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_pre_exp$delta_hct_10_0
## t = -22.249, df = 129, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -13.50271 -11.29729
## sample estimates:
## mean of x 
##     -12.4

Osmolality

osml_rmod_post_exp <- t.test(recovery_v_post_exp$delta_osml_10_8, 
                           mu = 0, alternative = "two.sided")
osml_rmod_post_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_post_exp$delta_osml_10_8
## t = 3.4782, df = 121, p-value = 0.0007021
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   4.330203 15.772256
## sample estimates:
## mean of x 
##  10.05123
osml_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_osml_10_0, 
                           mu = 0, alternative = "two.sided")
osml_rmod_pre_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_pre_exp$delta_osml_10_0
## t = 3.75, df = 130, p-value = 0.0002651
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   5.758213 18.618378
## sample estimates:
## mean of x 
##   12.1883

Group Export

Now, I put all the recovery t-test statistics into one dataframe, and export it.

recovery_stats <- broom.mixed::tidy(osml_rmod_pre_exp) %>%
  rbind(broom.mixed::tidy(osml_rmod_post_exp)) %>%
  rbind(broom.mixed::tidy(hct_rmod_pre_exp)) %>%
  rbind(broom.mixed::tidy(hct_rmod_post_exp)) %>%
  rbind(broom.mixed::tidy(SMI_rmod_pre_exp)) %>%
  rbind(broom.mixed::tidy(SMI_rmod_post_exp)) %>%
  mutate(response = c(rep("Plasma Osmolality (mmol/kg)", 2),
                      rep("Hematocrit (%)", 2),
                      rep("Body Condition (g')", 2)),
         pre_post_exp = c(rep(c("pre", "post"), 3)))
write.csv(recovery_stats, 
          "./results_statistics/recovery_stats.csv")

Figures

Colors & Shapes

Create custom colors, shapes, and labels to make it easy to apply changes across all figures.

CH_color <- brewer.pal(4, "Spectral")[c(4)]
HH_color <- brewer.pal(4, "Spectral")[c(2)]
CD_color <- brewer.pal(4, "Spectral")[c(3)]
HD_color <- brewer.pal(4, "Spectral")[c(1)]
my_colors <- c(CH_color, HH_color, CD_color, HD_color)

CH_shp <- 15
HH_shp <- 19
CD_shp <- 22
HD_shp <- 21
CH_shp_box <- 22
HH_shp_box <- 21
my_shapes <- c(CH_shp, HH_shp, CD_shp, HD_shp)
my_shapes_box <- c(CH_shp_box, HH_shp_box, CD_shp, HD_shp)

my_labels <- c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")

Body Mass

Means Only F5

ggplot() + 
  #plot these first so they end up on the "bottom"
  geom_smooth(data = dat_no_rehab,
                  aes(x = day_n,
                y = mass_g, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means,
                  aes(x = day_n,
                y = mean_mass, 
                color = tmt,
                group = tmt,
                ymin = mean_mass-se_mass, 
                ymax = mean_mass+se_mass),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_mass, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(seq(8,12)),
                     labels = c(seq(8,12))) +
  xlab("Day") + 
  ylab("Body Mass (g)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 10.8, unit = "pt")
) -> mass_fig_min
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
mass_fig_min
## `geom_smooth()` using formula = 'y ~ x'

Ending Values F6

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = mass_g, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = mass_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = mass_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_y_continuous(limits = c(7,15),
                     breaks = c(seq(7,15, by = 2)),
                     labels = c(seq(7,15, by = 2))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  xlab("") + 
  
  annotate(geom = "text", x = 4, y = 14.4, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 15, label = "BC", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 14.6, label = "AC", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 14.8, label = "A", 
           size = 3) +
  
  ylab("Body Mass (g)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              3.4), "mm")
        ) -> mass_end_boxplot
mass_end_boxplot

mass_emmeans
##                    tmt    emmean        SE       df  lower.CL upper.CL
## 1 Cool Humid (0.6 kPa) 10.667647 0.2517552 46.72727 10.161103 11.17419
## 2  Hot Humid (1.1 kPa)  9.727273 0.2559428 47.91893  9.212643 10.24190
## 3   Cool Dry (2.5 kPa) 10.581818 0.2559914 47.91893 10.067090 11.09655
## 4    Hot Dry (3.8 kPa)  9.547059 0.2521465 46.16198  9.039562 10.05456
##        response
## 1 Body Mass (g)
## 2 Body Mass (g)
## 3 Body Mass (g)
## 4 Body Mass (g)
mass_pairwise
##                                     contrast    estimate        SE       df
## 1 Cool Humid (0.6 kPa) - Hot Humid (1.1 kPa)  0.94037433 0.3575715 126.2004
## 2  Cool Humid (0.6 kPa) - Cool Dry (2.5 kPa)  0.08582888 0.3585976 127.4556
## 3   Cool Humid (0.6 kPa) - Hot Dry (3.8 kPa)  1.12058824 0.3554223 126.8612
## 4   Hot Humid (1.1 kPa) - Cool Dry (2.5 kPa) -0.85454545 0.3615808 127.8215
## 5    Hot Humid (1.1 kPa) - Hot Dry (3.8 kPa)  0.18021390 0.3583741 127.2163
## 6     Cool Dry (2.5 kPa) - Hot Dry (3.8 kPa)  1.03475936 0.3577985 126.4621
##      t.ratio    p.value      response
## 1  2.6298915 0.04674770 Body Mass (g)
## 2  0.2393459 0.99515381 Body Mass (g)
## 3  3.1528361 0.01074007 Body Mass (g)
## 4 -2.3633597 0.08951418 Body Mass (g)
## 5  0.5028653 0.95826222 Body Mass (g)
## 6  2.8920171 0.02307102 Body Mass (g)

SMI

Ind + Means

ggplot() + 
  geom_line(data = dat,
            aes(x = day_n,
                y = SMI, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means,
            aes(x = day_n,
                y = mean_SMI, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_SMI, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  geom_vline(xintercept = 9,
             linetype = "dashed",
             color = "darkgrey") +
  theme_classic() + 
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("Day") + 
  ylab("Body Condition (g)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0.1, #top
                             0.1, #right
                             0.1, #bottom
                             0.41 #left
                             ), "cm")) -> SMI_fig
SMI_fig

Means Only F5

ggplot() + 
  #plot these first so they end up on the "bottom"
  geom_smooth(data = dat_no_rehab,
                  aes(x = day_n,
                y = SMI, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means,
                  aes(x = day_n,
                y = mean_SMI, 
                color = tmt,
                group = tmt,
                ymin = mean_SMI-se_SMI, 
                ymax = mean_SMI+se_SMI),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_SMI, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(seq(8,12)),
                     labels = c(seq(8,12))) +
  xlab("Day") + 
  ylab("Body Condition (g')") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 10.8, unit = "pt")
) -> SMI_fig_min
SMI_fig_min
## `geom_smooth()` using formula = 'y ~ x'

LM + SE

ggplot() + 
  geom_smooth(data = dat_no_rehab,
              aes(x = day_n,
                  y = mass_g, 
                  color = tmt,
                  group = tmt),
              formula = y ~ x, 
              method = "lm", 
              se = T, 
              size = 1.2, 
              alpha = 0.2) + 
  theme_classic() + 
  scale_x_continuous(limits = c(0,8),
                     breaks = c(0, 2, 4, 6, 8)) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("Acclimation Day") + 
  ylab("Lizard Mass (g)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none"
) -> SMI_lm_fig
SMI_lm_fig

Ending Values F6

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = SMI, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = SMI_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = SMI_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_y_continuous(limits = c(7,15),
                     breaks = c(seq(7,15, by = 2)),
                     labels = c(seq(7,15, by = 2))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  xlab("") + 
  
  annotate(geom = "text", x = 4, y = 12.7, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 14.2, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 14.8, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 14.6, label = "A", 
           size = 3) +
  
  ylab("Body Condition (g')") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              3.4), "mm")
        ) -> SMI_end_boxplot
SMI_end_boxplot

SMI_emmeans
##                    tmt    emmean        SE       df  lower.CL upper.CL
## 1 Cool Humid (0.6 kPa) 10.773306 0.3101046 6.511477 10.028724 11.51789
## 2  Hot Humid (1.1 kPa)  9.624568 0.3116093 6.634876  8.879428 10.36971
## 3   Cool Dry (2.5 kPa) 10.654259 0.3115933 6.634043  9.909138 11.39938
## 4    Hot Dry (3.8 kPa)  9.420686 0.3101705 6.516199  8.676064 10.16531
##              response
## 1 Body Condition (g')
## 2 Body Condition (g')
## 3 Body Condition (g')
## 4 Body Condition (g')
SMI_pairwise
##                                     contrast   estimate        SE       df
## 1 Cool Humid (0.6 kPa) - Hot Humid (1.1 kPa)  1.1487389 0.2384728 126.0137
## 2  Cool Humid (0.6 kPa) - Cool Dry (2.5 kPa)  0.1190474 0.2391374 126.1449
## 3   Cool Humid (0.6 kPa) - Hot Dry (3.8 kPa)  1.3526208 0.2370101 126.0772
## 4   Hot Humid (1.1 kPa) - Cool Dry (2.5 kPa) -1.0296915 0.2411556 126.1941
## 5    Hot Humid (1.1 kPa) - Hot Dry (3.8 kPa)  0.2038818 0.2390104 126.1219
## 6     Cool Dry (2.5 kPa) - Hot Dry (3.8 kPa)  1.2335733 0.2385858 126.0337
##      t.ratio      p.value            response
## 1  4.8170640 2.420102e-05 Body Condition (g')
## 2  0.4978201 9.594278e-01 Body Condition (g')
## 3  5.7070166 4.643724e-07 Body Condition (g')
## 4 -4.2698224 2.211875e-04 Body Condition (g')
## 5  0.8530250 8.288557e-01 Body Condition (g')
## 6  5.1703556 5.287879e-06 Body Condition (g')

Hct

Ind + Means

ggplot() +
  geom_line(data = dat[complete.cases(dat$hematocrit_percent),],
            aes(x = day_n,
                y = hematocrit_percent, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means[complete.cases(means$mean_hct),],
            aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  geom_vline(xintercept = 9,
             linetype = "dashed",
             color = "darkgrey") +
  theme_classic() + 
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("") + 
  ylab("Hematocrit (%)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0.1, #top
                             0.1, #right
                             0.1, #bottom
                             0.41 #left
                             ), "cm")
        ) -> hct_fig
hct_fig
## Warning: Removed 12 rows containing missing values (`geom_point()`).

Means Only F5

ggplot() +
  geom_smooth(data = dat_no_rehab[complete.cases(dat_no_rehab$hematocrit_percent),],
                  aes(x = day_n,
                y = hematocrit_percent, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means[complete.cases(means$mean_hct),],
                  aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                group = tmt,
                ymin = mean_hct-se_hct, 
                ymax = mean_hct+se_hct),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means[complete.cases(means$mean_hct),],
            aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            alpha = 1, 
            fill = "white",
            size = 3) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(25, 30, 35, 40),
                     labels = c(25, 30, 35, 40),) +
  xlab("") + 
  ylab("Hematocrit (%)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 0, l = 10.8, unit = "pt")
        ) -> hct_fig_min
hct_fig_min
## `geom_smooth()` using formula = 'y ~ x'

Ending Values F6

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = hematocrit_percent, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = hct_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = hct_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  xlab("") + 
  scale_y_continuous(limits = c(10,55),
                     breaks = c(seq(10,50, by = 10)),
                     labels = c(seq(10,50, by = 10))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  
  annotate(geom = "text", x = 4, y = 47.6, label = "AB", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 44.6, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 54.6, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 54.6, label = "A", 
           size = 3) +
  
  ylab("Hematocrit (%)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              2.24), "mm")
        ) -> hct_end_boxplot
hct_end_boxplot
## Warning: Removed 3 rows containing missing values (`geom_point()`).

hct_emmeans
##                    tmt   emmean       SE       df lower.CL upper.CL
## 1 Cool Humid (0.6 kPa) 30.33352 1.808392 5.860156 25.88283 34.78422
## 2  Hot Humid (1.1 kPa) 25.77326 1.830146 6.142261 21.32008 30.22645
## 3   Cool Dry (2.5 kPa) 29.70335 1.815476 5.950456 25.25206 34.15464
## 4    Hot Dry (3.8 kPa) 27.63092 1.815279 5.948653 23.17978 32.08205
##         response
## 1 Hematocrit (%)
## 2 Hematocrit (%)
## 3 Hematocrit (%)
## 4 Hematocrit (%)
hct_pairwise
##                                     contrast   estimate       SE       df
## 1 Cool Humid (0.6 kPa) - Hot Humid (1.1 kPa)  4.5602630 1.272328 123.0398
## 2  Cool Humid (0.6 kPa) - Cool Dry (2.5 kPa)  0.6301769 1.254531 123.1123
## 3   Cool Humid (0.6 kPa) - Hot Dry (3.8 kPa)  2.7026074 1.253006 123.0665
## 4   Hot Humid (1.1 kPa) - Cool Dry (2.5 kPa) -3.9300861 1.288296 123.2298
## 5    Hot Humid (1.1 kPa) - Hot Dry (3.8 kPa) -1.8576556 1.283706 123.1026
## 6     Cool Dry (2.5 kPa) - Hot Dry (3.8 kPa)  2.0724305 1.262089 123.0607
##      t.ratio     p.value       response
## 1  3.5841894 0.002703439 Hematocrit (%)
## 2  0.5023206 0.958383983 Hematocrit (%)
## 3  2.1568990 0.141276348 Hematocrit (%)
## 4 -3.0506087 0.014665591 Hematocrit (%)
## 5 -1.4471033 0.472669199 Hematocrit (%)
## 6  1.6420637 0.359065321 Hematocrit (%)

Osml

Ind + Means

ggplot() + 
  geom_line(data = dat[complete.cases(dat$osmolality_mmol_kg_mean),],
            aes(x = day_n,
                y = osmolality_mmol_kg_mean, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means[complete.cases(means$mean_osml),],
            aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  geom_vline(xintercept = 9,
             linetype = "dashed",
             color = "darkgrey") +
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
  ylim(300,450) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("") + 
  ylab("Plasma Osmolality (mmol/kg)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme_classic() +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0.6, #top
                             0.1, #right
                             0.1, #bottom
                             0.1 #left
                             ), "cm")
        ) -> osml_fig
osml_fig
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 12 rows containing missing values (`geom_point()`).

Means Only F5

ggplot() + 
  geom_smooth(data = dat_no_rehab[complete.cases(dat_no_rehab$osmolality_mmol_kg_mean),],
                  aes(x = day_n,
                y = osmolality_mmol_kg_mean, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means,
                  aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                group = tmt,
                ymin = mean_osml-se_osml, 
                ymax = mean_osml+se_osml),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(seq(320,400, by = 10)),
                     labels = c(seq(320,400, by = 10))) +
  xlab("") + 
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 0, l = 1, unit = "pt")
        ) -> osml_fig_min
osml_fig_min
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 12 rows containing missing values (`geom_point()`).

Stats! Check Pairwise Diffs ~ Time

Since Plasma osmolality has such a nonlinear trend, we need to test whether the elevated values in the middle of the experiment are significantly different than the values taken before and/or after.

# first make sub-dfs for each tmt group
HH <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Hu")
HD <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Dr")
CH <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool H")
CD <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool D")

# next do pairwise tests for osml on the diff exp days, for each tmt group
pair_HH <- TukeyHSD(aov(data = HH, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_HD <- TukeyHSD(aov(data = HD, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_CH <- TukeyHSD(aov(data = CH, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_CD <- TukeyHSD(aov(data = CD, osmolality_mmol_kg_mean ~ day_factor)) #nonsig

# put into a df and export
osml_pairwise_df <- as.data.frame(pair_HD[[1]]) %>%
  rbind(as.data.frame(pair_HH[[1]])) %>%
  rbind(as.data.frame(pair_CD[[1]])) %>%
  rbind(as.data.frame(pair_CH[[1]])) %>%
  mutate(day_diff = paste("day", substr(rownames(.), 1, 3)),
         tmt = c(rep("Hot Dry",3),
                 rep("Hot Humid",3),
                 rep("Cool Dry",3),
                 rep("Cool Humid",3)),
         CI_95 = paste(round(lwr, digits = 2), round(upr, digits = 2), sep = ", "),
         diff = round(diff, digits = 2)) %>%
  dplyr::select(tmt, day_diff, diff, CI_95, p_adj = "p adj")

# save
write.csv(osml_pairwise_df, "./results_statistics/osmolality_pairwise_diffs.csv")

Nope, none of the differences between days within tmt groups are significantly different.

LM + SE

ggplot() + 
  stat_smooth(data = dat_no_rehab,
              aes(x = day_n,
                  y = osmolality_mmol_kg_mean, 
                  color = tmt,
                  fill = tmt,
                  group = tmt),
              formula = y ~ x, 
              method = "lm", 
              se = T, 
              size = 2, 
              alpha = 0.1) + 
  theme_classic() + 
  scale_x_continuous(limits = c(0,8), 
                     breaks = c(0, 2, 4, 6, 8)) +
  scale_color_brewer(palette = "Set2", name = "") +
  scale_fill_brewer(palette = "Set2", name = "") +
  xlab("Acclimation Day") + 
  ylab("Plasma Osmolality\n(mmol/kg)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 30),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 20),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none"
        ) -> osml_lm_fig
osml_lm_fig
## Warning: Removed 413 rows containing non-finite values (`stat_smooth()`).

Ending Values F6

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = osmolality_mmol_kg_mean, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = osml_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = osml_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_y_continuous(limits = c(290,470),
                     breaks = c(seq(300,450, by = 50)),
                     labels = c(seq(300,450, by = 50))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  xlab("") + 
  
  annotate(geom = "text", x = 4, y = 427, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 452, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 437, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 470, label = "A", 
           size = 3) +
  
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              0), "mm")
        ) -> osml_end_boxplot
osml_end_boxplot
## Warning: Removed 10 rows containing missing values (`geom_point()`).

osml_emmeans
##                    tmt   emmean       SE       df lower.CL upper.CL
## 1 Cool Humid (0.6 kPa) 350.3202 10.16423 4.874923 323.9893 376.6511
## 2  Hot Humid (1.1 kPa) 354.0963 10.13690 4.824373 327.7506 380.4420
## 3   Cool Dry (2.5 kPa) 349.1611 10.17830 4.903510 322.8413 375.4808
## 4    Hot Dry (3.8 kPa) 362.0244 10.15727 4.863220 335.6918 388.3569
##                      response
## 1 Plasma Osmolality (mmol/kg)
## 2 Plasma Osmolality (mmol/kg)
## 3 Plasma Osmolality (mmol/kg)
## 4 Plasma Osmolality (mmol/kg)
osml_pairwise
##                                     contrast   estimate       SE       df
## 1 Cool Humid (0.6 kPa) - Hot Humid (1.1 kPa)  -3.776107 5.031118 116.0115
## 2  Cool Humid (0.6 kPa) - Cool Dry (2.5 kPa)   1.159115 5.128594 116.0589
## 3   Cool Humid (0.6 kPa) - Hot Dry (3.8 kPa) -11.704175 5.078569 116.0347
## 4   Hot Humid (1.1 kPa) - Cool Dry (2.5 kPa)   4.935222 5.085161 116.0412
## 5    Hot Humid (1.1 kPa) - Hot Dry (3.8 kPa)  -7.928068 5.035558 116.0197
## 6     Cool Dry (2.5 kPa) - Hot Dry (3.8 kPa) -12.863291 5.114607 116.0115
##      t.ratio    p.value                    response
## 1 -0.7505502 0.87626200 Plasma Osmolality (mmol/kg)
## 2  0.2260104 0.99590771 Plasma Osmolality (mmol/kg)
## 3 -2.3046207 0.10288296 Plasma Osmolality (mmol/kg)
## 4  0.9705144 0.76645671 Plasma Osmolality (mmol/kg)
## 5 -1.5744172 0.39720576 Plasma Osmolality (mmol/kg)
## 6 -2.5150105 0.06282759 Plasma Osmolality (mmol/kg)

CEWL

Ind + Means

ggplot() + 
  geom_line(data = dat[complete.cases(dat$CEWL_g_m2h_mean),],
            aes(x = day_n,
                y = CEWL_g_m2h_mean, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  theme_classic() + 
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("Day") + 
  ylim(0,60) +
  ylab(bquote('CEWL (g/'*m^2*'h)')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "bottom"
        ) -> CEWL_fig
CEWL_fig

Means Only F3

ggplot() + 
  geom_errorbar(data = means[complete.cases(means$mean_CEWL),],
                  aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                group = tmt,
                ymin = mean_CEWL-se_CEWL, 
                ymax = mean_CEWL+se_CEWL),
                width = .1,
                #position=position_dodge(.1),
                alpha = 0.8) +
  geom_line(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                #linetype = tmt,
                group = tmt),
            alpha = 0.7, 
            size = .5) +
  geom_point(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(20, 25, 30, 35, 40),
                     limits = c(18,40)) +
  xlab("Day") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none"
        #legend.position = c(0.25,0.85)
        ) -> CEWL_fig_min
CEWL_fig_min

# use ggarrange so legend is centered
CEWL_fig_formatted <- ggarrange(CEWL_fig_min,
                                ncol = 1, nrow = 1,
                                common.legend = TRUE,
                                legend = "bottom")
# save
ggsave(filename = "experiment_CEWL_fig.pdf",
       plot = CEWL_fig_formatted,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 90)

LM + SE

ggplot() + 
  stat_smooth(data = dat,
              aes(x = day_n,
                  y = CEWL_g_m2h_mean, 
                  color = tmt,
                  fill = tmt,
                  group = tmt),
              formula = y ~ x, 
              method = "lm", 
              se = T, 
              size = 2, 
              alpha = 0.2) + 
  theme_classic() + 
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_color_brewer(palette = "Set2", name = "") +
  scale_fill_brewer(palette = "Set2", name = "") +
  xlab("Acclimation Day") + 
  ylab(bquote('CEWL (g/'*m^2*'h)')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 30),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 20),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none"
        ) -> CEWL_lm_fig
CEWL_lm_fig
## Warning: Removed 669 rows containing non-finite values (`stat_smooth()`).

delta CEWL ~ VPD F4

ggplot(data = dat_no_rehab_deltaCEWL) +
  geom_point(aes(x = VPD_kPa_tmttrial,
                   y = delta_CEWL, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 2,
              alpha = 0.4
             ) +
  geom_smooth(aes(x = VPD_kPa_tmttrial,
                   y = delta_CEWL), 
              se = F,
              formula = y ~ x,
              method = "lm",
              color = "black") +
  #geom_smooth(aes(x = VPD_kPa_tmttrial,
   #                y = delta_CEWL), 
    #          se = F,
     #         formula = y ~ poly(x, 2),
      #        method = "lm",
       #       color = "red") +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_color_manual(values = my_colors, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_y_continuous(limits = c(-13,40),
                     breaks = c(seq(-10,40, by = 10)),
                     labels = c(seq(-10,40, by = 10))
                     ) +
  scale_x_continuous(limits = c(0,4),
                     breaks = c(0.6, 1.1, 2.5, 3.8),
                     labels = c("0.6\nCH", 
                                "1.1\nHH", 
                                "2.5\nCD", 
                                "3.8\nHD")) +
  xlab("Vapor Pressure Deficit (kPa)") + 
  ylab(expression(Delta ~ 'CEWL')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "bottom",
        legend.spacing.y = unit(0, "mm")
        #plot.margin = unit(c(0, #top
         #                     0, #right
          #                    0, #bottom
           #                   0), "mm")
        ) -> CEWL_VPD_fig
CEWL_VPD_fig
## Warning: Removed 11 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 11 rows containing missing values (`geom_point()`).

# use ggarrange so legend is centered
CEWL_VPD_fig_formatted <- ggarrange(CEWL_VPD_fig,
                                ncol = 1, nrow = 1,
                                common.legend = TRUE,
                                legend = "bottom")
## Warning: Removed 11 rows containing non-finite values (`stat_smooth()`).
## Removed 11 rows containing missing values (`geom_point()`).
## Warning: Removed 11 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 11 rows containing missing values (`geom_point()`).
# save
ggsave(filename = "exp_CEWL_delta_VPD.pdf",
       plot = CEWL_VPD_fig_formatted,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 90)

Ending Values F6

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = CEWL_g_m2h_mean, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = CEWL_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = CEWL_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  scale_y_continuous(limits = c(9,63),
                     breaks = c(seq(10,60, by = 10)),
                     labels = c(seq(10,60, by = 10))) +
  xlab("") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) + 
  
  annotate(geom = "text", x = 4, y = 50, label = "C", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 63, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 41, label = "C", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 52, label = "A", 
           size = 3) +
  
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              0), "mm")
        ) -> CEWL_end_boxplot
CEWL_end_boxplot
## Warning: Removed 1 rows containing missing values (`geom_point()`).

CEWL_emmeans
##                    tmt   emmean       SE       df lower.CL upper.CL
## 1 Cool Humid (0.6 kPa) 30.36157 1.445712 10.68500 27.16810 33.55505
## 2  Hot Humid (1.1 kPa) 36.76428 1.446643 10.69577 33.56916 39.95941
## 3   Cool Dry (2.5 kPa) 23.72439 1.446548 10.69551 20.52946 26.91931
## 4    Hot Dry (3.8 kPa) 24.61233 1.434910 10.37491 21.43074 27.79392
##       response
## 1 CEWL (g/m2h)
## 2 CEWL (g/m2h)
## 3 CEWL (g/m2h)
## 4 CEWL (g/m2h)
CEWL_pairwise
##                                     contrast   estimate       SE       df
## 1 Cool Humid (0.6 kPa) - Hot Humid (1.1 kPa) -6.4027092 1.477747 125.0660
## 2  Cool Humid (0.6 kPa) - Cool Dry (2.5 kPa)  6.6371867 1.482776 125.4155
## 3   Cool Humid (0.6 kPa) - Hot Dry (3.8 kPa)  5.7492451 1.468760 125.2015
## 4   Hot Humid (1.1 kPa) - Cool Dry (2.5 kPa) 13.0398959 1.482912 125.4316
## 5    Hot Humid (1.1 kPa) - Hot Dry (3.8 kPa) 12.1519543 1.469693 125.2717
## 6     Cool Dry (2.5 kPa) - Hot Dry (3.8 kPa) -0.8879416 1.467098 125.0811
##      t.ratio      p.value     response
## 1 -4.3327496 1.740823e-04 CEWL (g/m2h)
## 2  4.4761900 9.849699e-05 CEWL (g/m2h)
## 3  3.9143543 8.436015e-04 CEWL (g/m2h)
## 4  8.7934385 9.026113e-14 CEWL (g/m2h)
## 5  8.2683648 1.057265e-12 CEWL (g/m2h)
## 6 -0.6052366 9.302438e-01 CEWL (g/m2h)

Exp CEWL ~ Osml

end_vals_CEWL_osml <- dat %>%
  dplyr::filter(day_n == 8) %>%
  dplyr::filter(complete.cases(CEWL_g_m2h_mean, osmolality_mmol_kg_mean))

ggplot(end_vals_CEWL_osml) + 
  aes(x = osmolality_mmol_kg_mean,
      y = CEWL_g_m2h_mean,
      color = tmt,
      shape = tmt) +
  geom_point(size = 1.5, 
             alpha = 0.8) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1, 
              alpha = 0.9) + 
  theme_classic() + 
  xlab(bquote('Osmolality (mmol '*kg^-1*')')) + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) + 
  scale_shape_manual(values = c(21,19, 22,15), name = "") +
  scale_fill_brewer(palette = "Spectral", name = "") +
  scale_color_brewer(palette = "Spectral", name = "") +
  scale_x_continuous(breaks = c(300, 350, 400, 450)) +
  scale_y_continuous(breaks = c(20, 30, 40, 50),
                     limits = c(12,57)) +
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.position = "bottom"
        ) -> exp_end_CEWL_osml_fig
exp_end_CEWL_osml_fig

ggsave(filename = "exp_CEWL_osml_fig.pdf",
       plot = exp_end_CEWL_osml_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 90)

Multi-Figures

over time:

ggarrange(osml_fig_min, 
          hct_fig_min, 
          mass_fig_min,
          ncol = 1, nrow = 3,
          labels = c("A", "B", "C"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> experiment_multi_fig
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 12 rows containing missing values (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 12 rows containing missing values (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#experiment_multi_fig

# export figure
ggsave(filename = "experiment_multi_fig.pdf",
       plot = experiment_multi_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 210)

end values:

ggarrange(CEWL_end_boxplot, 
          osml_end_boxplot, 
          hct_end_boxplot,
          mass_end_boxplot, 
          ncol = 2, nrow = 2,
          labels = c("A", "B", "C", "D"),
          widths = c(2, 2.045), heights = c(2, 2),
          common.legend = FALSE
          ) -> ending_values_multi_fig
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Warning: Removed 10 rows containing missing values (`geom_point()`).
## Warning: Removed 3 rows containing missing values (`geom_point()`).
#ending_values_multi_fig

ggsave(filename = "exp_end_val_multi_fig.pdf",
       plot = ending_values_multi_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 180, height = 150)